第三十六讲 R-多元线性回归中的交互作用及更优模型选择

在“R与生物统计专题”中，我们会从介绍R的基本知识展开到生物统计原理及其在R中的实现。以从浅入深，层层递进的形式在投必得医学公众号更新。

在前一讲中，我们介绍了如何用多个预测变量（x）建立用于预测连续型数值的结果变量（y）的多元线性回归模型（第三十五讲 R-多元线性回归）。

例如，要预测血糖，根据在胰岛素水平和怀孕次数，模型方程式为Glucose = b0 + b1*Insulin + b2*Pregnancies，其中b0是截距；b1 和b2是分别与预测变量血压、怀孕次数和怀孕次数相关联的回归系数。

上面的方程，也称为加法模型，仅研究预测变量的主要影响。它的前提条件是，假设给定的预测变量与结果之间的关系独立于其他预测变量。

但是，各个预测变量之间是不是相互独立的呢？如果怀孕次数的增加可以提高胰岛素水平对血糖的影响，则认为怀孕次数与胰岛素水平之间有协同效应，在统计学上，称为交互作用。

1. 多元线性回归中的交互作用方程式

具有两个预测变量（x1和x2）之间的相互作用的多元线性回归方程可写为：

y = b0 + b1*x1 + b2*x2 + b3*(x1*x2)

考虑我们的示例，它变为：

Glucose = b0 + b1*Insulin + b2* Pregnancies + b3*(Insulin* Pregnancies)

也可以写成：

Glucose = b0 + (b1 + b3* Pregnancies)*Insulin + b2* Pregnancies

或为：

Glucose = b0 + b1*Insulin + (b2 +b3*Insulin)* Pregnancies

b3 可以解释为胰岛素提高，而怀孕次数额外增加一个单位（反之亦然）。

2. R-多元线性回归中的交互作用

2.1 加载所需的R包

• tidyverse 用于数据可视化

• caret 用于简化机器学习的工作流程

library(tidyverse)
library(caret)

2.2 数据举例

我们将使用一组糖尿病的数据,它包含768人糖尿病相关的数据，包括怀孕情况，血糖，血压，皮肤厚度，胰岛素水平，怀孕次数，糖尿病谱系功能，年龄和糖尿病诊断结果（Outcome）。

注意

如需获取数据diabetes.csv，请关注投必得医学公众号，后台回复“diabetes.csv”获取数据。

导入数据：

my_data<-read.csv('diabetes.csv')

检查数据：

dim(my_data)

[1] 768   9

head(my_data)

输出结果

Pregnancies Glucose BloodPressure SkinThickness Insulin  PREGNANCIES DiabetesPedigreeFunction Age Outcome
1           6     148            72            35       0 33.6                    0.627  50       1
2           1      85            66            29       0 26.6                    0.351  31       0
3           8     183            64             0       0 23.3                    0.672  32       1
4           1      89            66            23      94 28.1                    0.167  21       0
5           0     137            40            35     168 43.1                    2.288  33       1
6           5     116            74             0       0 25.6                    0.201  30       0

数据清理

new_data<-my_data[my_data$Glucose>0 & my_data$Insulin >0,]
dim(new_data)

#[1] 393   9

研究问题：胰岛素水平和怀孕次数、及其交互作用预测血糖水平。

# 将数据分成训练数集和检验数集

set.seed(123)
training.samples <- new_data$Glucose %>%
createDataPartition(p = 0.8, list = FALSE)
train.data  <- new_data[training.samples, ]
test.data <- new_data[-training.samples, ]

2.3 加法模型

标准线性回归模型可以计算如下：

model1 <- lm(Glucose ~ Insulin + Pregnancies, data = train.data)
summary(model1)

输出结果

Call:
lm(formula = Glucose ~ Insulin + Pregnancies, data = train.data)


Residuals:
Min      1Q  Median      3Q     Max
-59.368 -17.110  -3.895  13.805  83.981


Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 93.91665    2.61539  35.909  < 2e-16 ***
Insulin      0.15017    0.01182  12.700  < 2e-16 ***
Pregnancies  1.55894    0.42649   3.655 0.000301 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 24.45 on 312 degrees of freedom
Multiple R-squared:  0.3731,  Adjusted R-squared:  0.369
F-statistic: 92.83 on 2 and 312 DF,  p-value: < 2.2e-16

# 用该模型计算验证数据集中预测值

predictions <- model1 %>% predict(test.data)

# 模型在验证数据集中的预测评估情况

# (a) 预测误差, RMSE
RMSE(predictions, test.data$Glucose)

[1] 25.86233

# (b) R平方
R2(predictions, test.data$Glucose)

[1] 0.3167198

2.4 交互效应

在R中，使用*运算符号来表示变量之间的交互作用：

model2 <- lm(Glucose ~ Insulin + Pregnancies + BloodPressure: Pregnancies,
data = new_data)

# 也可以使用以下R代码，两者输出一样

model2 <- lm(Glucose ~ Insulin* Pregnancies, data = train.data)
summary(model2)

输出结果

Call:
lm(formula = Glucose ~ Insulin * Pregnancies, data = train.data)


Residuals:
Min      1Q  Median      3Q     Max
-56.176 -16.922  -4.239  13.203  85.592


Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)         92.133076   3.301245  27.909   <2e-16 ***
Insulin              0.162520   0.018282   8.890   <2e-16 ***
Pregnancies          2.142610   0.784902   2.730   0.0067 **
Insulin:Pregnancies -0.003746   0.004229  -0.886   0.3763
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 24.45 on 311 degrees of freedom
Multiple R-squared:  0.3746,  Adjusted R-squared:  0.3686
F-statistic: 62.11 on 3 and 311 DF,  p-value: < 2.2e-16

# 用该模型计算验证数据集中预测值

# 模型在验证数据集中的预测评估情况

predictions <- model2 %>% predict(test.data)
RMSE(predictions, test.data$Glucose)

[1] 25.61056

# 计算 R平方

R2(predictions, test.data$Glucose)

[1] 0.3299708

2.5 结果解释

可以看出，当我们考虑胰岛素水平和怀孕次数的交互作用时，该交互作用并不能预测血糖水平，其系数与0间没有差异。即最终模型中，我们不应该考虑胰岛素水平和怀孕次数间的交互作用，其对血糖没有交互作用。

我们可以进一步使用anova()函数，对两个预测模型进行对比比较其R平方。

anova(model1,model2)

输出结果

Analysis of Variance Table
Model 1: Glucose ~ Insulin + Pregnancies
Model 2: Glucose ~ Insulin * Pregnancies
Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    312 186454
2    311 185984  1    469.35 0.7848 0.3763

从结果我们可以看出，model2并没有比model1好。残差平方和（Residual sum of squares，RSS）在model1中为186454，在model2中为185984。两者进行卡方检验，P = 0.2763，即两者残差平方和没有统计学差异。

所以，最后更好的模型为model1。

请注意

有时交互作用项很重要，但不是主要影响。如果模型中交互作用没有统计学意义，最终的模型中，不应该加入交互作用项。

如果模型中包括交互作用，并且交换作用有统计学意义，那么在模型中，即便原预测因素的系数与0没有统计学差异，也应该包括在最终的模型中。

参考内容：

1. Alboukadel Kassambara, Machine Learning Essentials: Practical Guide in R

好了，本期讲解就先到这里。小伙伴们赶紧试起来吧。

在之后的更新中，我们会进一步为您介绍R的入门，以及常用生物统计方法和R实现。欢迎关注，投必得医学手把手带您走入R和生物统计的世界。

提前预告一下，下一讲我们将学习多元线性回归中的多重共线性和方差膨胀因子。

当然啦，R语言的掌握是在长期训练中慢慢积累的。一个人学习太累，不妨加入“R与统计交流群”，和数百位硕博一起学习。

快扫二维码撩客服，

带你进入投必得医学交流群，

让我们共同进步！

↓↓

- END -

长按二维码关注「投必得医学」，更多科研干货在等你！